#!/usr/bin/perl -w
use strict;
use FindBin;
use Getopt::Long;
use lib "/net/cpp-group/Leo/bin";
use parse_fasta;
use Cwd;
use run_cmd;
use die_error;

#888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<"USAGE";

USAGE:

  codeml_parallel.pl 

      DESCRIPTION:
              Runs site specific codeml analyses in parallel on the cluster. 
  
      REQUIRED PARAMETERS:
              --cdna CDNA.FILE       N.B. These are file names not paths.
              --prot PEPTIDE.FILE          Set --working_dir if necessary
  
      OPTIONAL PARAMETERS:
              --f61                  Use full codon table instead of F3x4
                                     N.B. This is more biologically accurate 
                                          but often over-parameterises.
              --job_name  JOBNAME    Specify this to be able to 
                                          view job progress / kill jobs later
              --remove_gaps          Remove all columns which contain ANY
                                          missing residues (i.e. a gap)
              --msf                  Use alternative formats for peptide
              --clustal                   multiple alignment 
                                          (defaults to FASTA)
              --regex "(.*)"         Regular expression to find 
                                          corresponding peptide and cDNA
                                          identifiers. (See following.)
              --verbose              Print more information
              --working_dir              Print more information
              --help                 This screen
  

NOTES:
    Run this programme on the cluster head node (currently fgu205).
        (i.e. run "ssh fgu205" to logon etc.)
  

    The corresponding cDNA and peptide identifiers must match.
        By default, identifiers are the first letters until a space/tab.
        Otherwise, the identifiers can be selected using a regular 
        expression (in --regex)

CODEML
    Nested models 1 and 2 is a simplistic model.
    Nested models 7 and 8 uses a beta distribution and 
        may be over parameterised.
    Both should predict +ve selection to be trusted.

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    Make sure positively selected sites have a omega >> 1
       or they may merely be instances of neutral evolution.

    Always check for errors in the pipeline (All files ending in '.err')
       ***ESPECIALLY***  mapping.err
	   

RESULTS:
    See results.parameters

FAILURES:
    If it fails, it may be because the tree is too large/small.
        Check the pairwise dS values in 'codeml.pairwise/pairwise.dS'.
        The corresponding fitch-predicted dS tree is in 'fitch.tree'.

   Happy Hunting.
                  .... Leo (Goodstadt) 1/11/2004 to 17/08/2007

USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $remove_gaps = undef;
my $regex = '(\\S+)';

my $job_name = '';
my $working_dir;
my $prot_file;
my $cdna_file;
my $msf;
my $clustal;
my $f61;
my $verbose;
my $vverbose;
{
    open CMDLINE, ">cmd.line" or die "Error:\n\tCould not open cmd.line\n";
    print CMDLINE join (" ", @ARGV), "\n";
    close CMDLINE;
}
GetOptions(
			'prot=s'	=> \$prot_file,
			'cdna=s'	=> \$cdna_file,
			'regex=s'	=> \$regex,
			'msf'		=> \$msf,
			'clustal'	=> \$clustal,
			'f61'		=> \$f61,
			'verbose'	=> \$verbose,
			'vverbose'	=> \$vverbose,
			'help'		=> \$help,
			'remove_gaps'=> \$remove_gaps,
			'job_name=s'=> \$job_name,
			'working_dir=s'=> \$working_dir,
			);

die $usage if ($help);
die $usage unless ($prot_file && $cdna_file);

my $bin_dir = $FindBin::Bin;
$working_dir ||= getcwd();

die "Error\n\t:[$working_dir/$prot_file] does not exist.\n" unless (-e $prot_file);
die "Error\n\t:[$working_dir/$cdna_file] does not exist.\n" unless (-e $cdna_file);
my $verbose_cmd = $vverbose ? "--verbose" :"";
unless ($regex =~ /^".*"$/)
{
	$regex = '"'.$regex.'"';
}

$f61 = defined ($f61) ? 0 : 3;

$verbose = 1 if ($vverbose);














#_________________________________________________________________________________________
# 
# 	Safe make directory
# 
#_________________________________________________________________________________________
sub my_mkdir($)
{
	my ($dir_name) = @_;
	if (!-e $dir_name)
	{
		mkdir $dir_name
				or die "Error\n\tCould not create the directory [$dir_name]\n$!\n";
	}
	elsif (! -d  "$dir_name")
	{
		die "Error\n\t[$dir_name] is a file not a directory.\n$!\n";
	}

}


sub print_pairwise_ctl;

print "\n\n" if $verbose;




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#		Main logic


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#_________________________________________________________________________________________
# 
# 	Prepare sequences
# 
#_________________________________________________________________________________________
#
#	Prepare FASTA format protein sequence file
#
if ($msf)
{
	print
#STDERR
 "MSF --> FASTA ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/msf_to_fasta.pl $working_dir/$prot_file
		> $working_dir/$prot_file.fa
CMD
}
elsif ($clustal)
{
	print
#STDERR
 "CLUSTAL --> FASTA ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/clustal_to_fasta.pl $working_dir/$prot_file
		> $working_dir/$prot_file.fa
CMD
}
else
{
run_cmd(<<"CMD");
	cp $working_dir/$prot_file $working_dir/$prot_file.fa
CMD
}




#_________________________________________________________________________________________
# 
# 	Map cDNA onto aa
# 		
# 		Can run locally
# 
#_________________________________________________________________________________________
{
	print
#STDERR
 "Using protein alignment to align corresponding cDNA ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/map_cdna_onto_aa
		--peptide_sequences $working_dir/$prot_file.fa
        --cdna_sequences $working_dir/$cdna_file
        --peptide_output $working_dir/$prot_file.mapped.fa
        --cdna_output $working_dir/$cdna_file.mapped.fa
        --regex $regex
		--missing_identifiers $working_dir/missing.cdna_identifiers.err
        --err_log $working_dir/mapping.err
		--x_for_aa_mismatches X
		$verbose_cmd
CMD
	unlink ("$prot_file.fa");
}

if (-f "$working_dir/mapping.err")
{
	open_or_die(*ERROR, "$working_dir/mapping.err");
	print $_ for (<ERROR>);
}
if (-f "$working_dir/missing.cdna_identifiers.err")
{
	open_or_die(*ERROR, "$working_dir/missing.cdna_identifiers.err");
	print $_ for (<ERROR>);
}

#
#	Prepare cDNA onto aa
#
{
	print
#STDERR
 "Convert mapped cDNA to PAML format ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/seq_extract
		--sequences $working_dir/$cdna_file.mapped.fa
		--wrap_line 0
		--output_format p
		$verbose_cmd
		--e $working_dir/extract.err
		> $working_dir/seq.tmp
CMD
}

if (-f "$working_dir/extract.err")
{
	open_or_die(*ERROR, "$working_dir/extract.err");
	print $_ for (<ERROR>);
}

#_________________________________________________________________________________________
# 
#
#  -------Create pairwise Ks directory
#
#
#  -------make directories
#
my @start_val = ("0.1", "0.3", "1.3");
my_mkdir("codeml.pairwise");
for my $start_val(@start_val)
{
	print "Make directories for codeml analyses ...\n" if $verbose;
	my_mkdir ("$working_dir/codeml.$start_val");
}






#
#  -------Create PAML control files ..
#
print "Create PAML control files ...\n" if $verbose;
# pre-declaration
sub print_site_specific_ctl($);
for my $start_val(@start_val)
{
	print_site_specific_ctl($start_val);
}
print_pairwise_ctl();




#_________________________________________________________________________________________

#_________________________________________________________________________________________
#
#  -------Calculate pairwise Ks to build tree
#


{
	print "Calculate pairwise Ks to build tree ...\n" if $verbose;
	run_queue_cmd(<<"CMD", $job_name."_seq_pairs");
	$bin_dir/seq_pairs_kaks
		-aA
		-t $working_dir/codeml.pairwise
		-c $working_dir/$cdna_file.mapped.fa
		--err_log $working_dir/kaks.pairwise.err
		-p $working_dir/paml.err
		$verbose_cmd
		> $working_dir/kaks.results
CMD

if (-f "$working_dir/kaks.pairwise.err")
{
	open_or_die(*ERROR, "$working_dir/kaks.pairwise.err");
	print $_ for (<ERROR>);
}




print "Filter out Ks for build tree ...\n" if $verbose;
run_cmd(<<"CMD");
	$bin_dir/seq_pairs_filter_kaks_results
		-S
		-R $working_dir/filter.summary
		-r 0
		-f 0
		-s 0
		-t 0
		-n 0
		-l 0
		-e $working_dir/filter.err
		$verbose_cmd
		< $working_dir/kaks.results
		> $working_dir/codeml.pairwise/pairwise.dS
CMD

if (-f "$working_dir/filter.err")
{
	open_or_die(*ERROR, "$working_dir/filter.err");
	print $_ for (<ERROR>);
}


}



#
#  -------Build fitch tree from ks values
#
{
	print "Build fitch tree from ks values ...\n" if $verbose;

	# run fitch
	run_queue_cmd(<<"CMD", $job_name."_tree_fitch");
	$bin_dir/tree_fitch -i p
			$verbose_cmd
			< $working_dir/codeml.pairwise/pairwise.dS
			> $working_dir/fitch.tree
CMD
}



#
#  -------site specific analyses
#
{
	my @run_cmds;
	for my $i (0, 1, 2)
	{
		my $description = $verbose ? "Site specific analyses with a starting omega of $start_val[$i] ...\n": undef;
		my $dir = "$working_dir/codeml.$start_val[$i]";
	
		my $codeml_path = "$bin_dir/../../tools/paml/codeml ";
		my $ctl_file = "site_specific_kaks.$start_val[$i].ctl ";
		my $output_file = "> $working_dir/output.$start_val[$i]";
	
		my $exec_cmd = $codeml_path . $ctl_file . $output_file;
		push(@run_cmds,  [$exec_cmd, $job_name . "_" . $start_val[$i], $dir, $description]);
	}

	run_queue_parallel(\@run_cmds);
}


{
	print
#STDERR
 "Analyse results ...\n" if $verbose;

	run_cmd(<<"CMD");
	$bin_dir/parse_codeml_results.pl
			--para $working_dir/results.parameters
			--positive $working_dir/results.positive.fa
			$verbose_cmd
			$working_dir/site_specific.results.*
CMD
}


print "All finished...\n\n";
print "Results in results.parameters\n\n";
print "Nested models 1 and 2 is a simplistic model.\n";
print "Nested models 7 and 8 uses a beta distribution and may be over parameterised.\n";
print "I normally require all 3 to predict +ve selection before I believe anything.\n";
print "If it fails, it may be because the tree is too large/small.\n";
print "   Check the pairwise dS values in 'codeml.pairwise/pairwise.dS'.\n";
print "   The corresponding fitch-predicted dS tree is in 'fitch.tree'.\n";
print "Make sure any identified +vely selected sites have a large omega >> 1\n";
print "   or they may merely be instances of neutral evolution.\n";
print "Always check for errors in the pipeline (All files ending in '.err')\n";
print "   ***ESPECIALLY*** \n";
print "   mapping.err\n\n";
print "Happy Hunting.\n\n";
print "          .... Leo (Goodstadt) 27/05/2005\n\n\n";








#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	parameter files
#







sub print_pairwise_ctl
{
	open PAIR_CTL, ">$working_dir/codeml.pairwise/pairwise_kaks.ctl"
			or die "Error\n\t:Could not write to pairwise_kaks.ctl\n$!\n";
	print PAIR_CTL <<"CMD";
seqfile = ../seq.tmp
outfile = results
icode = 0	  			# universal code
noisy = 0
verbose = 0
runmode = -2  			# pairwise
seqtype = 1	  			# codons
CodonFreq = $f61		# F3X4; Could be 3 for codon table Was error fixed thanks to Andreas
aaDist = 0	  			# ??? equal Could be geometric or linear or G1974 Miyata
model = 0	  			# one dN/dS over tree
NSsites = 0	  			# one dN/dS over sequence
Mgene = 0		    	# ??? rates; 1: separate
fix_omega = 0
omega = .4
fix_alpha = 1		 	# gamma distribution shape set to 1 for single rate
alpha = 0		     	# gamma distribution shape set to 1 for single rate
fix_rho = 1
rho = 0
Malpha = 0		     	# ??? different alphas for genes
ncatG = 1		     	# # of categories in dG of NSsites models5D
clock = 0		     	# ??? no clock 1:clock
getSE = 1
RateAncestor = 1	 	# ??? rates alpha > 0 or ancestral rates 1 or 2
Small_Diff = .5e-6
method = 0		     	# simultaneous, 1: one branch at a time
fix_blength = 0		 	# ignore branch lengths from supplied ree
fix_kappa = 0
kappa = 2
CMD
	if ($remove_gaps)
	{
		print PAIR_CTL <<"CMD";
cleandata = 1
CMD
	}
	else
	{
		print PAIR_CTL <<"CMD";
cleandata = 0
CMD
	}

}

sub print_site_specific_ctl($)
{
	my ($start_omega) = @_;
	my $file = ">$working_dir/codeml.$start_omega/site_specific_kaks.$start_omega.ctl";
	open FITCH_CTL, $file
			or die "Error\n\t:Could not write to $file\n$!\n";
print FITCH_CTL <<"CMD";
seqfile = ../seq.tmp
treefile = ../fitch.tree
outfile = ../site_specific.results.$start_omega
icode = 0				* universal code
noisy = 0
verbose = 1
runmode = 0				* user supplied tree
seqtype = 1				* codons
CodonFreq = $f61		* 2 = F3x4,	3 = F64
aaDist = 0				* ??? equal Could be geometric or linear or G1974 Miyata
model = 0				* one dN/dS over tree
NSsites = 1 2 7 8
						* 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
						* 5:gamma;6:2gamma;7:beta;8:beta&w;9:beta&gamma;
						* 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
						* 13:3normal>0
						* one dN/dS over sequence
Mgene = 0				* ??? rates; 1: separate
fix_omega = 0
omega = $start_omega
fix_alpha = 1			* gamma distribution shape set to 1 for single rate
alpha = 0				* gamma distribution shape set to 1 for single rate
fix_rho = 1
rho = 0
Malpha = 0				* ??? different alphas for genes
ncatG = 10				* * of categories in dG of NSsites models5D
clock = 0				* ??? no clock 1:clock
getSE = 1
RateAncestor = 0		* ??? rates alpha > 0 or ancestral rates 1 or 2
Small_Diff = .5e-6
method = 0				* simultaneous, 1: one branch at a time
fix_blength = 0			* ignore branch lengths from supplied ree
fix_kappa = 0
kappa = 2CMD
CMD
	if ($remove_gaps)
	{
		print FITCH_CTL <<"CMD";
cleandata = 1
CMD
	}
	else
	{
		print FITCH_CTL <<"CMD";
cleandata = 0
CMD
	}
}


sub print_fitch_parameters
{
	my $file = ">$working_dir/fitch.parameters";
	open FITCH_CTL, $file
			or die "Error\n\t:Could not write to $file\n$!\n";
	print FITCH_CTL <<"CMD";
build_tree.dS
0
7
L
J
33
50
2
Y
CMD
	close FITCH_CTL;
}

